

%This function calculates country average markups...
function [markup, profits] = avgMarkups(m, rhohat)

markup = zeros(4,1);

%Find break between germany and denmark in rho vector
ger_rhos_end = (m.num_firms_ger+1)*m.num_pro_ger;

% For active firms, the (additive) markup is (1/(1-rho)) we wish to calculate
%  the expected markup which will weight each firm by its win probability.
%  However, we must correct for the fringe, which does not price
%  strategically.
weighted_markup = rhohat ./ (1 - rhohat);

%  Here we create a mask to "zero out" the fringe markup...the fringe is the first firm in every project
fringe_mask = [ kron(ones(m.num_pro_ger,1), [0; ones(m.num_firms_ger,1)]); ... German projects
                kron(ones(m.num_pro_dnk,1), [0; ones(m.num_firms_dnk,1)]) ];   %Danish projects
%apply the mask...
weighted_markup = weighted_markup .* fringe_mask;

%These are the expected country average markups paid equally weighted
%across projects
markup(1) =  sum(weighted_markup(1:ger_rhos_end))/m.num_pro_ger;
markup(2) =  sum(weighted_markup(ger_rhos_end+1:end))/m.num_pro_dnk;

%These are the expected country average markups paid weighted by megawatt output.
markup(3) = sum(kron(m.ger_mw,ones(m.num_firms_ger+1,1)).*weighted_markup(1:ger_rhos_end))/m.num_pro_ger;
markup(4) = sum(kron(m.dnk_mw,ones(m.num_firms_dnk+1,1)).*weighted_markup(ger_rhos_end+1:end))/m.num_pro_dnk;

%% Generate estimates of profits for each firm in each market

p_matrix_ger    = reshape(weighted_markup(1:(m.num_firms_ger+1)*m.num_pro_ger)',m.num_firms_ger+1,m.num_pro_ger)' ;
p_matrix_dnk    = reshape(weighted_markup((m.num_firms_ger+1)*m.num_pro_ger+1:end)',m.num_firms_dnk+1,m.num_pro_dnk)' ;

profit_matrix_ger = repmat(m.ger_mw,1,m.num_firms_ger).*p_matrix_ger(:,2:end);
profit_matrix_dnk = repmat(m.dnk_mw,1,m.num_firms_dnk).*p_matrix_dnk(:,2:end);

%We assume that fewer firms are active in denmark and zero-pad the
%difference here....this will work in the "usual case" and for the "Add 1"
%counterfactuals...however it is NOT general.
profits = [ sum(profit_matrix_ger); ...
            sum(profit_matrix_dnk) zeros(1, m.num_firms_ger - m.num_firms_dnk)];


